Gas-Liquid Coexistence in the Primitive Model for Water 
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We evaluate the location of the gas-liquid coexistence line and of the associated critical point 
for the primitive model for water (PMW), introduced by Kolafa and Nezbeda [J. Kolafa and I. 
Nezbeda, Mol. Phys. 61, 161 (1987)]. Besides being a simple model for a molecular network 
forming liquid, the PMW is representative of patchy proteins and novel colloidal particles interacting 
with localized directional short-range attractions. We show that the gas-liquid phase separation 
is metastable, i.e. it takes place in the region of the phase diagram where the crystal phase is 
thermodynamically favored, as in the case of particles interacting via short-range attractive spherical 
potentials. Differently from spherical potentials, we do not observe crystallization close to the critical 
point. The region of gas-liquid instability of this patchy model is significantly reduced as compared 
to equivalent models of spherically interacting particles, confirming the possibility of observing 
kinetic arrest in an homogeneous sample driven by bonding as opposed to packing. 

PACS numbers: 61.20. Ja,61.20.Gy,61.20.Ne 



I. INTRODUCTION 

This article presents a detailed numerical study of the 
critical point and gas-liquid coexistence of a model intro- 
duced several years ago by Kolafa and Nezbeda [H as a 
primitive model for water (PMW). The water molecule 
is described as a hard-sphere with four interaction sites, 
arranged on a tetrahedral geometry, which are meant to 
mimic the two hydrogen protons and the two oxygen lone- 
pairs of the water molecule. The PMW has been studied 
in detail in the past, since it is both a valid candidate 
for testing theories of bond association 0, H, 0, H, H, 0, Hi 
as well as a model able to re pro duce the thermodynamic 
anomalies of water [ij, |g, llfj, llll, LL2J ■ The Wertheim the- 
ory @) 9 has been carefully compared to numerical stud- 
ies, suggesting a good agreement between theoretical pre- 
dictions and numerical data in the temperature T region 
where no significant ring formation is observed [ill [l2| . 
More recently, the slow dynamics of this model has been 
studied using a newly developed code for event-driven dy- 
namics of patchy particles [121 ]. It has been shown that, 
at low density, there are indications of a gas-liquid phase 
separation. At intermediate density the system can be 
cooled down to the smallest temperatures at which equi- 
libration is feasible with present day computational facil- 
ities without any signature of phase separation. In this 
region, dynamics progressively slows down following an 
Arrhenius law, consistently with the expected behavior 
of strong network- forming liquids. 

This primitive model, besides its interest as an elemen- 
tary model for water, is also representative of a larger 
class of models which are nowadays receiving consider- 
able interest, namely models of particles interacting with 
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localized directional short-range attractions. This class 
of models includes, apart from network forming molec- 
ular systems, also proteins Hp . [Til [l5l liH ] and newly 
designed colloidal particles [17]] • Recently, the phase di- 
agram of patchy particles has been studied [181 ] in the 
attempt to estimate the role of the surface patchiness in 
dictating their dynamic and thermodynamic behavior. It 
has been suggested that the number of possible bonds M, 
as opposed to the fraction of surface with attractive inter- 
actions, is the key ingredient in determining the width of 
the unstable region of the phase diagram [1^, [2(| Hl| . A 
study of the evolution of the critical point on decreasing 
M shows a clear monotonic trend [181 ] and, more impor- 
tantly, in the direction of decreasing critical packing frac- 
tion. Such a trend is not observed in spherical potentials 
on decreasing the attraction range. Thus, in low-valence 
particle systems, a window of packing fraction (f> values 
opens up in which it is possible to reach very low T (and 
hence states with extremely long bond lifetimes) without 
encountering phase separation. This favors the establish- 
ment of a spanning network of long-living bonds, which 
in the colloidal community provides indication of gel for- 
mation but which, in the field of network forming liquids, 
would be rather classified as glass formation [221 ]. 

As previously mentioned, we present here an accu- 
rate estimate of the location of the critical point for the 
PMW (for which M — 4), based on finite-size scaling, 
and of the associated gas-liquid coexistence phase dia- 
gram. Comparing with the known fluid-crystal phase 
coexistence loci [llj we are able to show that the gas- 
liquid phase separation is metastable as compared to the 
crystal state. Despite its metastability, we have never ob- 
served crystallization, not even close to the critical point, 
an observation which could be of relevance in the case 
of the crystallization of patchy proteins or colloids. We 
also show that the number density p of the liquid phase 
at coexistence is comparable to the crystal (diamond) 
density, much smaller than the characteristic liquid den- 
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FIG. 1: Pictorial representation of the model studied. Each 
particle is modeled as an hard-core sphere (grey large sphere 
of diameter a). The four interaction sites are located in a 
tetrahedral arrangement. Two of the sites (the H-sites, in 
green) are located on the surface whereas the remaining two 
sites (LP, in brown) are located inside the sphere, at distance 
0.45cr from the center. Only H-sites and L-sites on different 
particles interact with a square well interaction of depth uo 
and width 8 — 0.15cr. 

sity observed in systems of particles interacting through 
spherically symmetric potentials. 

II. THE MODEL 

In the PMW, each particle is composed by a hard 
sphere of diameter a and four additional sites arranged 
according to a tetrahedral geometry. Two of these sites, 
the proton sites H, are located on the surface of the hard 
sphere, i.e. at distance 0.5a from the center of the parti- 
cle, while the two remaining sites, the lone-pair sites LP, 
are placed within the hard sphere at a distance 0.45cr 
from its center. Besides the hard-sphere interaction, pre- 
venting different particles to sample relative distances 
smaller than <r, only the H and LP sites of distinct par- 
ticles interact via a square well (SW) potential usw of 
width 6 = 0.15a and depth uq, i.e. 

usw(r) = | r > 5 (1) 

where r is here the distance between H and LP sites. 
We remark that the choice S = 0.15<r guarantees that 
multiple bonding can not take place at the same site. 

III. SIMULATIONS 

Three types of Monte Carlo simulations (Metropo- 
lis Grand Canonical, GCMC, Umbrella Sampling Grand 
Canonical, US-GCMC and Gibbs Ensemble GEMC) have 
been performed, with the aim of locating the critical 



point of the model, assessing the consistency with the 
Ising universality class and evaluating the gas-liquid co- 
existence curve. Throughout the remaining sections, we 
use uo as energy scale, a as length scale and reduced 
units in which fee = 1, thus measuring T and the chemi- 
cal potential p, in units of uo/ks and u$ respectively. 

We define a MC step as an average of N/\ attempts to 
translate and rotate a randomly chosen particle and an 
average Nn attempts to insert or remove a particle. The 
translation in each direction is uniformly chosen between 
±0.05 and the rotations are performed around a random 
axis of an angle uniformly distributed between ±0.1 rad. 
Unless otherwise stated, Na/Nn = 500. The choice of 
such a large ratio between translation/rotation and in- 
sertion/deletion attempts is dictated by the necessity of 
ensuring a proper equilibration at fixed N. This is partic- 
ularly important in the case of particles with short-range 
and highly directional interactions, since the probability 
of inserting a particle with the correct orientation and 
position for bonding is significantly reduced as compared 
to the case of spherical interactions. 



A. Critical Point Estimate and Finite Size Scaling 
Analysis 

The location of the critical point has been performed 
through the comparison of the fluctuation distribution of 
the ordering operator M. at the critical point with the 
universal distribution characterizing the Ising class [23| . 
The ordering operator M. of the gas-liquid transition is 
a linear combination A4 ~ p + su, where p is the number 
density, u is the energy density of the system, and s is 
the field mixing parameter. Exactly at the critical point, 
fluctuations of M. are found to follow a known universal 
distribution, i.e. apart from a scaling factor, the same 
that characterizes the fluctuations of the magnetization 
in the Ising model (23|. Finite Size Scaling (FSS) anal- 
ysis has been performed to test the appartenance of the 
PMW model to the three dimensional Ising class. Recent 
applications of this method to soft matter can be found 
in Ref. [Ullilli]. 

To locate the critical point we perform, at each size, 
GCMC simulations, i.e. at fixed T, /i and volume V. We 
follow the standard procedure of tuning T and /x until 
the simulated system shows ample density fluctuations, 
signaling the proximity of the critical point. Once a rea- 
sonable guess of the critical point in the (T, /x)-plane has 
been reached, we start at least 10 independent GCMC 
simulations to improve the statistics of the fluctuations 
in the number of particles ./V in the box and of the po- 
tential energy E. 

The precise evaluation of the critical temperature and 
chemical potential at all system sizes is performed by 
a fitting procedure associated to histogram reweighting. 
We briefly recall [27[ that this technique allows us to pre- 
dict, from the joint distribution P(N, E;T, fi) obtained 
from a simulation, the distribution at a different T" and 
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// through the ratio of the Boltzmann factors as follows: 



e -fiE pPfiN 



e (P-P')E e (P'n'-l3»,)N _ 



P(N,E;n',T) 

(2) 

where the proportionality constant can be calculated im- 
posing the normalization of P(N, E; \J ,T'). The re- 
weighting procedure offers reliable results provided T" 
and // are within few percents of T and /i. 

We implement a least-square fit procedure to evaluate 
the values of T, /i and s for which the reweighted distri- 
bution of M. is closest to the known form for Ising-like 
systems [28|. The result of the fit provides the best esti- 
mate for the critical temperature T c , the critical chemical 
potential /i c and s. The location of the critical point has 
been performed for boxes of sizes 6 through 9, corre- 
sponding to an average number of particles ranging from 
about 60 to over 400. The simulations for larger boxes 
have been started at the critical parameters calculated 
for the closer smaller box. The L = 9 simulation re- 
quired about one month of CPU time for each of the 10 
studied replicas on a 2.8 GHz Intel Xeon. 



B. Phase coexistence 

We have calculated the phase coexistence curve us- 
ing US-GCMC and GEMC simulations. The US-GCMC 
method is a natural extension of the standard GCMC 
method used to locate the critical point discussed above. 
Once the density fluctuations at the critical point have 
been evaluated, one can again apply the histogram 
reweigthing method to predict the shape of the density 
fluctuations at T < T c for any fixed n value. The value 
of [i for which coexistence between a gas and a liquid 
phase is present is chosen by selecting the value /z x for 
which the areas underneath the two peaks of P(M), the 
ordering operator fluctuation distribution, are equal. 

We stress that performing a standard (Metropolis) 
GCMC simulation at a temperature even a few percent 
lower than T c is not feasible, due to the large free-energy 
barrier separating the two phases which would prevent 
the simulated system to sample both liquid and gas con- 
figurations, thus yielding physically meaningless results. 
Several techniques have been developed to overcome this 
problem, among which the Umbrella Sampling Monte 
Carlo [29| which has been used to perform the calcu- 
lations presented herein. The US is a biased sampling 
technique which aims to flatten the free energy barrier 
between the two phases modifying the standard GCMC 
insertion/removal [30( probabilities as follows: 



jGCMC w(N) 



ins ~~ ins w(N+l) 



(3) 



P l 



= P 



rem w(N— 1) 



In these last equations w(N) is a forecast on the real 
P(N) which we have repeatedly extracted from previous 
higher T simulations through the histogram reweighting 



technique. If the predicted fluctuation distribution w(N) 
is a good approximation to the real P(N), the resulting 
biased fluctuation will result flat in N and the system will 
thus not experience any difficulty in crossing the barrier 
between the liquid and gas phase. As shown in [29j, the 
unbiased fluctuation distribution can be easily recovered 
adding a reverse bias to the results obtained with the 
insertion/removal probability in Eq. [3] 

Starting from the critical point, we evaluate the phase 
diagram iterating the above procedure, progressively low- 
ering T, down to the point where equilibration was not 
achieved any longer. Indeed, on cooling, due to the for- 
mation of a well-connected tetrahedral network, the dy- 
namics in the liquid side slows down considerably 
We have been very careful in progressively increasing the 
ratio N^/Nn to compensate for the slowing down of the 
dynamics and the extremely slow equilibration times of 
the liquid phase. At the lower T, N A /N N = 10000. This 
sets a bound to the smallest T which can be investigated. 

We have also implemented a Gibbs Ensemble evalu- 
ation of the phase diagram. The GEMC method was 
designed 31| to study coexistence in the region where 
the gas-liquid free-energy barrier is sufficiently high to 
avoid crossing between the two phases. Since nowadays 
this is a standard method in computational physics, we 
do not discuss it here and limit ourselves to noting that 
also in this case it is important to progressively increase, 
on cooling, the ratio N&/Nn to account for the slow dy- 
namics characterizing the liquid state. We have studied a 
system of (total) 350 particles which partition themselves 
into two boxes whose total volume is 2868<r 3 , correspond- 
ing to an average density of p — 0.122. At the lowest T 
this corresponds to roughly 320 particles in the liquid 
box (of side w 8a) and about 30 particles in the gas box 
(of side w 13c). Equilibration at the lowest reported T 
required about three months of computer time. 



IV. RESULTS 

We start showing the distributions of the ordering 
operator fluctuations at the (apparent) critical points 
for all studied sizes. Implementing the fitting proce- 
dure described in Sec. IIII A[ using histogram reweight- 
ing, we first evaluate the values of the critical param- 
eters and the shape of the density fluctuations at the 
critical point. The resulting distributions, for all inves- 
tigated L, are shown as a function of the scaled variable 
x = a~^L$l v {M - M c ) in Fig. \21 Here v is the critical 
exponent of the correlation length and j3 is critical ex- 
ponent of the order parameter. Within the d = 3 Ising 
universality class v — 0.629 and (3 — 0.326. From the fits 
we find that the non-universal amplitude is = 0.34 
(independent on L). The size dependence of T c and fi c 
for L = 6 to L = 9 is shown in Fig. [3] Finite size scal- 
ing predicts T c ~ L-C«+i)/" and p c ~ i-C^ 1 )/", where 
9 = 0.54 (Ising) is the universal correction to the scal- 
ing exponent [23|]. Fig. [3Ja-b) shows that the size de- 
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FIG. 2: Distributions of the fluctuations of the ordering op- 
erator for all studied sizes at the apparent critical point. The 
factor aj^ = 0.34 has been chosen to scale P(x) to unit vari- 
ance. The theoretical curve for the Ising model (full line) is 
from 12811 . 



pendence of the critical parameters is consistent with 
the expected universality class. Extrapolating the ob- 
served size dependence to L — > oo it is possible to pro- 
vide an estimate for the bulk behavior of the PMW po- 
tential. We find T bulk = 0.1083(3), p b c ulk = -1.265(1). 
Fig-EI c ) shows the L dependence of p*, the latter being 
the system density (evaluated via histogram reweighting 
for each size) at T bulk and p, bulk . Finite size scaling pre- 
dicts p* ~ L-^- 1 /"). The resulting L dependence of p* 
is consistent with the theoretical prediction, despite the 
fact that p* data are the noisiest, since the estimate of 
density is the most delicate and the error is at least of the 
order of one particle over the volume. The field mixing 
parameter s has been difficult to infer precisely due to 
its small value and to the discrete nature of the square 
well interactions which makes the distribution of N and 
E quantized over integer numbers. A value of s m 0.07, 
consistent with the results reported in for a similar 
model, allows a simultaneous fit of all distribution func- 
tions, independently from L. 

Next we discuss the gas-liquid coexistence curve for 
the PMW. Since only close to the critical point size ef- 
fects are relevant, we have performed the calculations at 
L = 6. Fig. [4] shows the behavior of the density fluctua- 
tions P{p) (in log scale), evaluated with US-GCMC sim- 
ulations (see Sec. MI B]) along the coexistence curve. The 
difference between the peak and the valley in ln[P(p)] is 
a measure, in unit of k^T, of the activation free-energy 
Fbarrier I k-g,T needed to cross from the gas to the liquid 
phase and vice-versa. At the lowest studied temperature, 
the barrier reaches a value of about 20fceT, which would 
clearly be impossible to overcome without the use of a 
biased sampling technique. Fig. 2] also shows the T de- 
pendence of the barrier height, which indeed becomes of 
the order of the thermal energy close to the critical point. 




-(d-l/v) 

FIG. 3: Size dependence of the apparent critical temperature 
T C (L) (a), the apparent critical chemical potential p c (L) (b) 
and the density at the true critical point p*(L) (c). 



The phase diagram resulting from our calculations is 
reported in Fig. 0(a) . Both US-GCMC and GEMC data 
are reported, showing a perfect agreement for all stud- 
ied T proving that, despite the long equilibration times 
required, an accurate determination of the phase dia- 
gram for this model of patchy particles can nowadays be 
achieved. Fig. 0(b) shows the same data, together with 
the fluid-crystal coexistence lines calculated by Vega and 
Monson and complemented with the bond percola- 
tion line from Ref. [12]. We also report in the graph a 
so-called iso-diffusivity line, i.e. the set of points in the 
T — p plane in which the diffusion coefficient D is con- 
stant. By selecting the smallest value of D which can be 
calculated with the present time computational facilities, 
the iso-diffusivity line provides an estimate of the shape 
of the glass line. Fig. [5] also shows the coexistence curve 
evaluated from the Wertheim theory, from Ref. fill ] . 

Several observations arise from the data in Fig. [5] First 
of all, we observe that the critical point, as well as the 
liquid branch of the coexistence curve, are characterized 
by a percolating (but transient) structure of bonds be- 
tween pairs of H and LP. This is not unexpected, since 
the propagation of infinite range correlations, character- 
istic of a critical point, does require the presence of a 
spanning cluster [32|]. Particle diffusion is still signifi- 
cant, despite the presence of the transient percolating 
network of bonds, as shown by the comparison between 
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FIG. 4: (a): Logarithm of the distribution of density fluctu- 
ations P(p) at T = 0.1082, 0.1060, 0.1049, 0.1037, 0.1026, 
0.1009, 0.0999, providing the negative of the density depen- 
dence of the free energy. The difference between the value of 
P(p) at the valley and at the top is a measure of the free en- 
ergy barrier separating the gas and the liquid state (b): Free 
energy barrier (in units of feT) between gas and liquid versus 
T at coexistence. 
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FIG. 5: (a): Calculated gas- liquid coexistence for the PMW. 
Both GEMC (diamonds) and US-GCMC (circles) simulations 
results are shown. The star indicates the bulk (L — > co) 
critical point estimate, (b): Extended phase diagram for the 
PMW, including: (i) the theoretical (Wertheim) gas-liquid co- 
existence line (dashed), (ii) the field of stability of the fluid, of 
the high-density crystal (HDS) and of the low-density crys- 
tal (LDS) phases from Ref. [ll|; (iii) an iso-diffusivity line 
(for the smallest investigated value of D) from Ref. jl2l ]; (iv) 
percolation line from Ref. [12J. 



the phase-coexistence and the small D iso-diffusivity line. 
In the region of T where dynamics becomes so slow that 
equilibration can not be achieved on the present compu- 
tational time scale, it becomes impossible also to evaluate 
the gas-liquid coexistence. 

The gas-liquid coexistence is found to be metastable 
respect to fluid-crystal coexistence, in analogy with the 
case of particles interacting through spherical short-range 
potentials. This notwithstanding and differently from the 
case of spherical interacting particles, we never observe 
crystallization close to the critical point. This suggests 
that the increase in local density brought in by the critical 
fluctuations does not sufficiently couple with the orien- 
tational ordering required for the formation of the open 
diamond crystal structure. 

The comparison between the GEMC and US-GCMC 
simulation phase diagram with the theoretical predic- 
tions of the Wertheim theory suggests that the latter 
provides a quite accurate estimate for T c , whereas the p 



dependence is only approximate. Indeed, the Wertheim 
theory predicts a vapor-liquid critical point at T c — 
0.1031 and p c = 0.279 [11]. These differences between 
theoretical predictions and numerical data confirm the 
conclusions that have been previously reached for mod- 
els of patchy particles with different number of sticky 
spots [18[. 

Finally, we note that in a very limited T-interval (less 
than 10 per cent of T c ), the liquid density approxima- 
tively reaches a value comparable to the diamond crys- 
tal density, which for the present model has been calcu- 
lated [HI 0.4880 < p < 0.6495. Thus, as for the crystal 
state, the density of the liquid phase at coexistence is 
rather small, approximatively a factor of two smaller as 
compared to the case of spherically interacting particles. 

Finally, for the sake of completeness we show in Fig. [S] 
the location of the gas-liquid coexistence in the /z — T 
plane. 
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FIG. 6: Chemical potential (circles) and fugacity z = 
(diamonds) versus temperature T at coexistence for L = 6. 



V. CONCLUSIONS 

In this article we have presented an accurate determi- 
nation of the critical point and of the gas-liquid phase 
coexistence curve for a primitivejnodel for water, intro- 
duced by Kolafa and Nezbeda Despite its original 
motivation, the PMW can also be studied as a model for 
patchy colloidal particles and, perhaps, as an elementary 
model for describing patchy interactions in proteins. To 
this extent, it is particularly important to understand the 
qualitative features of the phase diagram, the stability or 
metastability of the gas-liquid line and the propensity to 
crystallize. 

We have shown that the critical fluctuations are consis- 
tent with the Ising universality class, both via the analy- 
sis of the shape of the density fluctuation distribution and 
via the size dependence of the critical parameters. The 
determination of the critical point allows us to prove that, 
for the present model, the liquid phase is not present in 
equilibrium, and it is only observed in a metastable con- 
dition. In the case of spherical interaction potentials, the 
thermodynamic stability of the gas-liquid critical point 
is achieved when the range of the interaction potential 
becomes of the order of one tenth of the particle diame- 



ter [33| . Data reported in this article (see Fig. [5]) confirm 
that the property of short-range potentials of missing a 
proper equilibrium liquid phase is retained in short-range 
patchy particles. It would be of interest to study for 
which critical value of the interaction range a proper liq- 
uid phase appears in such patchy particles. 

The evaluation of the gas-liquid coexistence has been 
performed using two distinct methods, the US-GCMC 
and the GEMC, and a very good agreement has been 
recorded despite the difficulties in equilibrating the liq- 
uid phase at low T. The present results, together with 
previous studies of the crystal phase and of the slow low- 
T dynamics, offer a coherent picture of the behavior of 
the model. It is shown that phase separation is interven- 
ing only at low T for p < 0.6, a value significantly smaller 
than the corresponding value for spherical interaction po- 
tentials. This is consistent with the fact that the number 
of attractive interactions in which a particle can be en- 
gaged (four in the present model) controls the minimum 
density of the liquid phase. These results elucidate the 
different nature of bonded liquids (the so-called network 
forming) with respect to the category of spherically in- 
teracting liquids. For bonded systems, there is an inter- 
mediate region, between the high density packed struc- 
ture and the unstable region, in which gas-liquid phase 
separation is observed, which is not accessible to spher- 
ically interacting particles. In this intermediate density 
region, the system is structured in a percolating network 
and both static and dynamic quantities are controlled by 
the presence of bonds. This picture, which has been de- 
veloping in a progression of studies [HI, 0j|, [22|, |34j, is 
confirmed by the present results. 
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